Large Seedling Linear Regression
Kathleen Stutzman
2024-01-10
Libraries:
# libraries -------------------
library(dplyr)
library(tidyverse)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(TMB)
library(emmeans)
library(DHARMa)
library(lmtest)
library(ggeffects)
library(performance)Set working directory
knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')Functions:
# functions -------------------Source large seedling data:
source("Scripts/LS_Import.R")Pivot wider to create dataframe where each row is for one plot and has total details for each species group
ls_merge2 <- ls_merge %>%
select(Plot_No, Region, Treat_Type, Site, Species_Groups, Total) #this is dropping browse and stump sprout data
ls_merge2 <- ls_merge2 %>%
pivot_wider(names_from = Species_Groups, values_from = Total)Import time since data and add it to the large seedling dataset
time_since <- read_csv("CleanData/Treat_Year_Data.csv")
ls_merge3 <- merge(ls_merge2, time_since, by = 'Site')
#log transform time from treatment data
ls_merge3$l.TFT <- log(ls_merge3$Time_from_Treat)Run the ‘Add_BA’ script and merge with dataset:
source("Scripts/Add_BA.R")
# merge with ls dataset -------------------
ls_merge4 <- merge(ls_merge3, prism_BA, by = "Plot_No")Run ‘Ground_Data.R’ script and add it to large seedling dataset:
source("Scripts/Ground_Data.R")
# merge with ls dataset -------------------
ls_merge5 <- merge(ls_merge4, ground3, by = "Plot_No")Source and add veg data
source("Scripts/Veg_Data.R")
# merge with ls dataset
ls.all <- merge(ls_merge5, veg3, by = "Plot_No")
rm(ls_merge5,
ls_merge2,
ls_merge3,
ls_merge4)The large seedling count data is taken in 10m2 plots; basal area is measured in hectares; veg and soil data is taken in 1m2 plots.
Large seedling data will be converted into 1m2 plots in order to compare across and reduce the amount of scales of data collection to two: 1m2 plots and per hectare observations.
Create log transformed categories for newly added variables, then select for just the desired variables:
ls.all$PIRI.1m <- ls.all$PIRI/10
ls.all$SO.1m <- ls.all$Shrub_Oak/10
ls.all$Other.1m <- ls.all$Other/10
ls.all$l.PIRI1 <- log(ls.all$PIRI.1m + 1)
ls.all$l.SO1 <- log(ls.all$SO.1m + 1)
ls.all$l.other1 <- log(ls.all$Other.1m + 1)
ls.all2 <- ls.all %>%
select(Treat_Type, Region, Site, Plot_No, PIRI, PIRI.1m, l.PIRI1, Shrub_Oak, SO.1m, l.SO1, Other, Other.1m, l.other1, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>%
arrange(Treat_Type)Select just for numerical vs log and then look at paired plots:
#not transformed
ls.num <- ls.all2 %>%
select(PIRI, Shrub_Oak, Other, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)
ggpairs(ls.num)
ggpairs(ls.num, aes(color = Treat_Type))
#log transformed
ls.numl <- ls.all2 %>%
select(l.PIRI, l.SO, l.other, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)
ggpairs(ls.numl)
ggpairs(ls.numl, aes(color = Treat_Type))
rm(ls.num,
ls.numl)
rm(ls.all2)Can see the correlation coefficients for linear (Pearsons) relationships. None of them appear very strong, except for ones that are analogs (avg LD vs mineral soil exposure; ba/ha vs piri ba/ha)
Log transformed average litter depth and basal area per hectare have a weak relationship (corr 0.33), which does make sense.
The above script a data naming conventions are the same as in the LS_GLMM script
PIRI Modeling
I tried with variables at 1m scale - didn’t work well. I also did model elimination not using treatment type - that also did not work well
Revised data set with LS observations at 10m2 scale; this means no non-interger values for the poisson distro
ls.all3 <- ls.all
ls.all3$l.PIRI <- log(ls.all3$PIRI + 1)
ls.all3$l.SO <- log(ls.all3$Shrub_Oak + 1)
ls.all3$l.other <- log(ls.all3$Other + 1)
ls.all3 <- ls.all3 %>%
select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI, Shrub_Oak, l.SO, Other, l.other, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>%
arrange(Treat_Type)testing zero inflation
tapply(ls.all3$PIRI, ls.all3$Region, summary)## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01639 0.00000 1.00000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.0137 0.0000 1.0000
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 0.0 0.0 0.5 0.0 14.0
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.266 0.000 9.000# no PIRI in ME
tapply(ls.all3$PIRI, ls.all3$Treat_Type, summary) ## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.04274 0.00000 3.00000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.09804 0.00000 4.00000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 1.304 1.250 14.000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.1103 0.0000 9.0000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.02299 0.00000 1.00000I started modeling with poisson distro and excluding treatment type. Models without treatment type failed in model fit. Going to run variable elimination again, starting with models with treatment type
ls.m13 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = nbinom1)
AIC(ls.m13) #288.9## [1] 288.9363ls.m14 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = nbinom1)
AIC(ls.m14) #287.8## [1] 287.8049ls.m15 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = nbinom1)
AIC(ls.m15) #285.9## [1] 285.9361ls.m16 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = nbinom1)
AIC(ls.m16) #284.4## [1] 284.4039ls.m17 <- glmmTMB(PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = nbinom1)
AIC(ls.m17) #283.9## [1] 283.9195lrtest(ls.m13, ls.m14, ls.m15, ls.m16, ls.m17) #non significant## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral +
## avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l +
## l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 3: PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 4: PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 5: PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 15 -129.47
## 2 14 -129.90 -1 0.8686 0.3513
## 3 13 -129.97 -1 0.1312 0.7172
## 4 12 -130.20 -1 0.4678 0.4940
## 5 11 -130.96 -1 1.5156 0.2183ls.m18 <- glmmTMB(PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = nbinom1)
AIC(ls.m18) #284.4## [1] 284.3545lrtest(ls.m17, ls.m18) #p = 0.12## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + l.SO + avgLD_l + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 |
## Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -130.96
## 2 10 -132.18 -1 2.4349 0.1187ls.m19 <- glmmTMB(PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = nbinom1)
summary(ls.m19) #282.4## Family: nbinom1 ( log )
## Formula:
## PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1 | Site/Plot_No)
## Data: ls.all3
##
## AIC BIC logLik deviance df.resid
## 282.4 320.3 -132.2 264.4 489
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 7.864e-09 8.868e-05
## Site (Intercept) 8.389e+00 2.896e+00
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Dispersion parameter for nbinom1 family (): 1.84
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.0011 1.8849 -4.245 2.19e-05 ***
## Treat_TypeFallRx 2.8946 1.8244 1.587 0.112610
## Treat_TypeHarvest 5.8075 1.9747 2.941 0.003271 **
## Treat_TypeMowRx 2.2997 1.7665 1.302 0.192968
## Treat_TypeSpringRx 2.9793 1.9726 1.510 0.130956
## avgLD_l -0.9013 0.2489 -3.622 0.000293 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lrtest(ls.m18, ls.m19) #0.85## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + avgLD_l + l.Veg_Total + offset(l.TFT) + (1 |
## Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 10 -132.18
## 2 9 -132.19 -1 0.0353 0.851# this works, but ls.m19 has lower AIC
ls.m20 <- glmmTMB(PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
ziformula = ~Region,
family = nbinom1)
AIC(ls.m20) #286.6## [1] 286.6293rm(ls.m13, ls.m14, ls.m15, ls.m16, ls.m17, ls.m18)Ok, test model fit with treat type (same as model 11)
#Best model
ls.m19 <- glmmTMB(PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = nbinom1)
summary(ls.m19) #282.4## Family: nbinom1 ( log )
## Formula:
## PIRI ~ Treat_Type + avgLD_l + offset(l.TFT) + (1 | Site/Plot_No)
## Data: ls.all3
##
## AIC BIC logLik deviance df.resid
## 282.4 320.3 -132.2 264.4 489
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 7.864e-09 8.868e-05
## Site (Intercept) 8.389e+00 2.896e+00
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Dispersion parameter for nbinom1 family (): 1.84
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.0011 1.8849 -4.245 2.19e-05 ***
## Treat_TypeFallRx 2.8946 1.8244 1.587 0.112610
## Treat_TypeHarvest 5.8075 1.9747 2.941 0.003271 **
## Treat_TypeMowRx 2.2997 1.7665 1.302 0.192968
## Treat_TypeSpringRx 2.9793 1.9726 1.510 0.130956
## avgLD_l -0.9013 0.2489 -3.622 0.000293 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ls.m19_sr <- simulateResiduals(ls.m19, n = 1000, plot = TRUE) #passes
testResiduals(ls.m19_sr) #pases
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.037938, p-value = 0.4705
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0031118, p-value = 0.814
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01711847
## sample estimates:
## outlier frequency (expected: 0.00118473895582329 )
## 0## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.037938, p-value = 0.4705
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.0031118, p-value = 0.814
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01711847
## sample estimates:
## outlier frequency (expected: 0.00118473895582329 )
## 0testZeroInflation(ls.m19_sr) #passes
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99996, p-value = 0.972
## alternative hypothesis: two.sidedtestQuantiles(ls.m19_sr) #passes
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 0.3178
## alternative hypothesis: both#This is a better fit (using nbinom instead of poisson). Convergence issues for nbinom2. Would different distribution change the variables included? - No, I went back and ran model elimination with nbinom1 distro instead of poisson. same variables are significant
ls.m20_sr <- simulateResiduals(ls.m20, n = 1000, plot = TRUE) #passes
testResiduals(ls.m20_sr) #passes
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.028424, p-value = 0.8159
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.13922, p-value = 0.774
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.56
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.006024096
## sample estimates:
## outlier frequency (expected: 0.000863453815261044 )
## 0.002008032## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.028424, p-value = 0.8159
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.13922, p-value = 0.774
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.56
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.006024096
## sample estimates:
## outlier frequency (expected: 0.000863453815261044 )
## 0.002008032testZeroInflation(ls.m20_sr) #passes
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.98999, p-value = 0.608
## alternative hypothesis: two.sidedtestQuantiles(ls.m20_sr) #passes
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 0.4708
## alternative hypothesis: bothTrying pairwise comparison of model fit:
emmeans(ls.m19, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.000518 0.000953 Inf 1.41e-05 0.0191
## FallRx 0.009362 0.016220 Inf 3.14e-04 0.2793
## Harvest 0.172374 0.280718 Inf 7.08e-03 4.1946
## MowRx 0.005164 0.008595 Inf 1.98e-04 0.1348
## SpringRx 0.010190 0.017281 Inf 3.67e-04 0.2829
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 0.0553 0.10093 Inf 1 -1.587 0.5060
## Control / Harvest 0.0030 0.00593 Inf 1 -2.941 0.0271
## Control / MowRx 0.1003 0.17716 Inf 1 -1.302 0.6901
## Control / SpringRx 0.0508 0.10026 Inf 1 -1.510 0.5557
## FallRx / Harvest 0.0543 0.10307 Inf 1 -1.535 0.5395
## FallRx / MowRx 1.8129 3.06288 Inf 1 0.352 0.9967
## FallRx / SpringRx 0.9187 1.74321 Inf 1 -0.045 1.0000
## Harvest / MowRx 33.3770 61.40768 Inf 1 1.907 0.3137
## Harvest / SpringRx 16.9152 33.56315 Inf 1 1.425 0.6112
## MowRx / SpringRx 0.5068 0.93049 Inf 1 -0.370 0.9960
##
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scaleControl vs Harvest is the only significant different (p = 0.0271)
Would things fit better with zi ~Region? Answer: NO, good model fit but higher AIC
GG Plot for ls piri
# PIRI LS -------------------
ls.p1 <- ggpredict(ls.m19, terms = c("avgLD_l", "Treat_Type"))
ggplot(ls.p1, aes(x, predicted, color = group))+
geom_line(linewidth = 1.25, show.legend = FALSE) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
theme(axis.text = element_text(size = 16), axis.title = element_text(size = 20))+
labs(x = "Average leaf litter depth \n (log transformed)",
y = "Pitch pine stems \n (adjusted for time)")
ggsave(filename = "LS_PIRI_avgld.tiff", path = "Plots/PIRI_GLMM", width = 7, height = 5, device = "tiff", dpi = 700)Shrub oak models
A losing proposition. No models pass fit tests.
Abandon all hope.
Are SO LS present in all regions and treatment types?
Yes
tapply(ls.all3$Shrub_Oak, ls.all3$Region, summary)## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 0.00 10.50 18.29 34.00 68.00
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 4.452 6.000 30.000
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 0.00 8.00 14.49 19.00 73.00
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 8.00 13.00 16.42 23.50 47.00
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 6.25 17.00 19.14 27.00 74.00tapply(ls.all3$Shrub_Oak, ls.all3$Treat_Type, summary)## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 2.000 5.684 9.000 47.000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 2.25 13.00 16.08 25.00 67.00
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.500 5.768 10.000 32.000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 13.00 26.00 27.68 40.00 74.00
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 0.0 7.0 12.6 17.0 61.0Tried a lot of distributions. Models are generally underdispersed. Have used poisson, nbinom, genpois, & compois distros. Using zi ~1 helps.
# so.ls1 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
# data = ls.all3,
# ziformula = ~1,
# family = genpois)
#AIC(so.ls1) # fails with genpois
# test piri ba
so.ls2 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
ziformula = ~1,
family = genpois)
AIC(so.ls2) #3266.7
#lrtest(so.ls1, so.ls2) #p = 0.3
# test total BA
so.ls3 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
ziformula = ~1,
family = genpois)
AIC(so.ls3) #3268
lrtest(so.ls2, so.ls3) # p = 0.06
# test mineral
so.ls4 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
ziformula = ~1,
family = genpois)
AIC(so.ls4) #3268.6
lrtest(so.ls3, so.ls4) # p = 0.1
# test avg ld
so.ls5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
ziformula = ~1,
family = genpois)
AIC(so.ls5) #3267.3
lrtest(so.ls5, so.ls4) # p = 0.4
# test piri*
so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
ziformula = ~1,
family = genpois)
AIC(so.ls6) #fails
lrtest(so.ls6, so.ls5) # p = 0.008**
# test other*
so.ls7 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.Veg_Total + l.PIRI + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
ziformula = ~1,
family = genpois)
AIC(so.ls7) #fails
lrtest(so.ls5, so.ls7) # p = 0.003
# test veg*
so.ls8 <- glmmTMB(Shrub_Oak ~ Treat_Type +l.other + l.PIRI + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
ziformula = ~1,
family = genpois)
AIC(so.ls8) # 3271.1
# test piri ba, model failed at start
so.ls9 <- glmmTMB(Shrub_Oak ~ Treat_Type +l.other + l.PIRI + l.Veg_Total + l.BA_piri + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
ziformula = ~1,
family = genpois)
AIC(so.ls9) #3268.4
lrtest(so.ls9, so.ls5) # p = 0.4Best Model
so.ls5 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
ziformula = ~1,
family = genpois)
summary(so.ls5)## Family: genpois ( log )
## Formula:
## Shrub_Oak ~ Treat_Type + l.PIRI + l.other + l.Veg_Total + offset(l.TFT) +
## (1 | Site/Plot_No)
## Zero inflation: ~1
## Data: ls.all3
##
## AIC BIC logLik deviance df.resid
## 3267.3 3317.8 -1621.6 3243.3 486
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 0.08187 0.2861
## Site (Intercept) 0.97499 0.9874
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Dispersion parameter for genpois family (): 8.01
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.94791 0.42709 -6.902 5.12e-12 ***
## Treat_TypeFallRx 3.20832 0.48348 6.636 3.22e-11 ***
## Treat_TypeHarvest 2.40848 0.55615 4.331 1.49e-05 ***
## Treat_TypeMowRx 4.32943 0.44516 9.726 < 2e-16 ***
## Treat_TypeSpringRx 3.65532 0.48789 7.492 6.78e-14 ***
## l.PIRI -0.38785 0.19308 -2.009 0.0446 *
## l.other -0.11699 0.05036 -2.323 0.0202 *
## l.Veg_Total 0.17110 0.07100 2.410 0.0160 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6709 0.3736 -7.149 8.73e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(so.ls5)## [1] 3267.271so.ls5_sr <- simulateResiduals(so.ls5, n = 1000, plot = T)
testResiduals(so.ls5_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.033648, p-value = 0.6257
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.10303, p-value = 0.002
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.009086345
## sample estimates:
## outlier frequency (expected: 0.00122489959839357 )
## 0## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.033648, p-value = 0.6257
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.10303, p-value = 0.002
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.009086345
## sample estimates:
## outlier frequency (expected: 0.00122489959839357 )
## 0testDispersion(so.ls5_sr, alternative = "less")
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.10303, p-value = 0.001
## alternative hypothesis: lesscheck_collinearity(so.ls5) #low## # Check for Multicollinearity
##
## * conditional component:
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
## Treat_Type 1.04 [1.00, 1.43] 1.02 0.96 [0.70, 1.00]
## l.PIRI 1.02 [1.00, 2.41] 1.01 0.98 [0.42, 1.00]
## l.other 1.02 [1.00, 3.06] 1.01 0.98 [0.33, 1.00]
## l.Veg_Total 1.03 [1.00, 1.68] 1.01 0.97 [0.60, 1.00]Pairwise with treatment type
emmeans(so.ls5, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')FallRx, Harvest, MowRx, & SpringRx all significantly different from Control. Only Harvest & MowRx significantly different from each outher
testing other distros
so.ls6 <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = poisson)
AIC(so.ls6)
so.ls6a <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = genpois)
AIC(so.ls6a)
# compois distro fails, won't produce DHARMa results
so.ls6b <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = tweedie)
AIC(so.ls6b)
so.ls6b_sr <- simulateResiduals(so.ls6b, n = 1000, plot = TRUE)
testDispersion(so.ls6b_sr) #fails
so.ls6c <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = t_family)
AIC(so.ls6c)
so.ls6c_sr <- simulateResiduals(so.ls6c, n = 1000, plot = TRUE)
testResiduals(so.ls6c_sr)
testDispersion(so.ls6c_sr) #passes
testQuantiles(so.ls6c_sr)
so.ls6d <- glmmTMB(Shrub_Oak ~ Treat_Type + l.other + l.BA_HA + l.Mineral + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
data = ls.all3,
family = gaussian)
AIC(so.ls6d)
so.ls6d_sr <- simulateResiduals(so.ls6d, n = 1000, plot = TRUE) #fails
plot(so.ls6d_sr)
testResiduals(so.ls6d_sr)hist(ls.all3$Shrub_Oak)
ggplot(ls.all3, aes(x=Shrub_Oak))+
geom_histogram(binwidth = 2)+
facet_grid(rows = vars(Treat_Type))
Graph model
# SO LS -------------------
ls.s1 <- ggpredict(so.ls5, terms = c("l.PIRI", "Treat_Type"))
ggplot(ls.s1, aes(x, predicted, color = group))+
geom_line(linewidth = 1.25) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
theme(axis.text = element_text(size = 16), axis.title = element_text(size = 20))+
labs(x = "Count of pitch pine seedlings (log)",
y = "Shrub oak stems \n(adjusted for time)")
#ggsave(filename = "LS_PIRI_avgld.tiff", path = "Plots/PIRI_GLMM", width = 7, height = 5, device = "tiff", dpi = 700)
ls.s2 <- ggpredict(so.ls5, terms = c("l.other", "Treat_Type"))
ggplot(ls.s2, aes(x, predicted, color = group))+
geom_line(linewidth = 1.25, show.legend = FALSE) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
theme(axis.text = element_text(size = 16), axis.title = element_text(size = 20))+
labs(x = "Count of other seedlings (log)",
y = "Shrub oak stems \n(adjusted for time)")
# -------------------
ls.s3 <- ggpredict(so.ls5, terms = c("l.Veg_Total", "Treat_Type"))
ggplot(ls.s3, aes(x, predicted, color = group))+
geom_line(linewidth = 1.25, show.legend = FALSE) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
theme(axis.text = element_text(size = 16), axis.title = element_text(size = 20))+
labs(x = "Total understory vegetative cover (log)",
y = "Shrub oak stems \n(adjusted for time)")
# low collinearity but still gives me a little pause..Interesting - it is always Mowing -> SpringRx -> Fall Rx -> Harvest -> Control
I’ve tried poisson, nbinom1, nbinom1 with ZI, nbinom2, genpois, compois, tweedie, gaussian and none of these distros work. They are all underdispersed
I’ve deleted more models than you’ve imagined.
Genpois with zi ~1 seems like the closest I’m going to get